# generate with matlab:
# for i = 1 : 20
# 	phi = randn(60);
#	eulpaths(i, 1:61) = InitrecursiveMonteCarlo(0.06, 1/30, 2, 0.2, phi, 100);
#   intpaths(i, 1:61) = InitintegralMonteCarlo(0.06, 1/30, 2, 0.2, phi, 100);
# end
# save backward_euler.res eulpaths -ASCII
# save integral.res intpaths -ASCII


paths <- read.table("backward_euler.res")
pathsint <- read.table("integral.res")
walks <- as.data.frame(t(paths))
plot(walks$V1,ylim=c(min(paths),max(paths)), type="n",xlab="Time", ylab="Stock value")
for (i in seq(1, length(walks))) {
	lines(walks[,i],col="blue",lty=9)
	}

title(paste(length(walks), "realizations of a Stock price using an Euler approximation and Integral approximation"))

par(new=T)

walksint <- as.data.frame(t(pathsint))
plot(walksint$V1,ylim=c(min(pathsint),max(pathsint)), type="n",xlab="",ylab="",axes=F)
for (i in seq(1, length(walksint))) {
	lines(walksint[,i],col="dark gray",lty=9)
	}
lines(as.data.frame(mean(pathsint[,])),col="black")
par(new=T)
plot(walks$V1,ylim=c(min(paths),max(paths)), type="n",xlab="",ylab="",axes=F)
lines(as.data.frame(mean(paths[,1:length(paths)])),col="blue")
legend(1,200,c("Euler approximation", "Integral approximation", "mean of the Euler approximations", "mean of the Integral approximations"), col=c("blue", "dark gray", "blue","black"), lty=c(9, 9, 1, 1))